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Abstract 

As fractional diffusion equations can describe the early breakthrough and 
the heavy-tail decay features observed in anomalous transport of contami- 
nants in groundwater and porous soil, they have been commonly employed 
in the related mathematical descriptions. These models usually involve long- 
time range computation, which is a critical obstacle for its application, im- 
provement of the computational efficiency is of great significance. In this 
paper, a semi-analytical method is presented for solving a class of time- 
fractional diffusion equations which overcomes the critical long-time range 
computation problem of time fractional differential equations. In the proce- 
dure, the spatial domain is discretized by the finite element method which re- 
duces the fractional diffusion equations into approximate fractional relaxation 
equations. As analytical solutions exist for the latter equations, the burden 
arising from long-time range computation can effectively be minimized. To 
illustrate its efficiency and simplicity, four examples are presented. In addi- 
tion, the method is employed to solve the time-fractional advection-diffusion 
equation characterizing the bromide transport process in a fractured granite 
aquifer. The prediction closely agrees with the experimental data and the 
heavy-tail decay of anomalous transport process is well-represented. 

Key words: Anomalous transport, Mittag-Leffier function, finite element 
method, time-fractional diffusion equation 
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1. Introduction 



For the contaminant transport processes in soil and groundwater, diffu- 
sion equations (such as diffusion equation, advection-dispersion equation and 
advection-reaction-diffusion equation) are the traditional governing equations 
[H, 0, 0, 0] • I n the past several decades, however, more and more evidences 
show that some of the critical features in contaminant transport through 
complex porous media cannot be described by the conventional diffusion 
equations 



12, 13 



These features include early breakthrough 
and heavy-tail decay of the contaminant as well as the scale-dependent co- 
efficients 



10, 11 



They have led to the increasing use of the fractional 
diffusion equations for modeling contaminant transport in porous media. 

The theoretical research on fractional diffusion equation models has re- 
ceived considerable success in physical modeling and experimental result 
analysis in the last decade 



15, 16 . Nevertheless, numerical meth- 



ods for solving fractional diffusion equations are still immature for practical 
applications in which the spatial problem domains are geometrically com- 
plex and large whilst the long-time range predictions are often desired. The 
main obstacle is that the fractional derivative has its global nature, compared 
with the locally expressed classic derivatives, the computational cost of time 
fractional derivative term will increase dramatically with time evolution. It 
brings computational challenge of approximating fractional order equations 
with the finite difference or the finite element methods 17J, [18j, |34( . Though 



short memory method is proposed to tackle the long-time range computa- 
tion, it has been proved to bring accuracy deterioration in many cases 0,0 



Moreover, when computing high Peclet number problems, numerical schemes 



may produce oscillating solutions [22j, |23|, |24], |25j, |26j. The problem is more 



annoying in fractional advection-dispersion equations. 

Until now, the finite difference method has been widely employed for 
solving anomalous diffusion equations with successes in short-time and small 
spatial scale problems 27, 28[. Since the finite element method (FEM) is 
more suitable for modeling large and geometrically complicated spatial do- 
mains, the investigation of the FEM for fractional diffusion type equations 
has attracted much attention in recent years. Roop et al. approximated the 
solutions of steady state space-fractional advection-dispersion equations in 



18, 29 



two spatial dimensions using the Galerkin and least-squares FEMs 
Huang et al. proposed an unconditionally stable FEM approach to solve 
the one-dimensional space-fractional advection-dispersion equation, and suc- 



2 



cessfully applied it to simulate the atrazine transport in a saturated soil 
column Deng developed a FEM for the numerical resolution of the 
space and time-fractional Fokker-Planck equation with its convergence order 



is 0(k 2 ~ a + h^), a and fi are time and spatial derivative orders 30j. Zhuang 
et al. investigated the Galerkin finite element approximation of symmetric 
space-fractional partial differential equations, and proved the stability and 



convergence of the proposed schemes [3J|. Zheng et al. discussed the FEM 



for the space-fractional advection diffusion equation with non-homogeneous 



initial-boundary condition [32l . l33l . l34j . Though all aforementioned works in- 



dicate that FEMs play an important and increasing role in the applications of 
fractional diffusion type equation models, the efficient and simple numerical 
methods for fractional diffusion equations are still urgently needed. 

In this paper, we introduce a semi-analytical FEM for time-fractional 
diffusion equations which can be expressed in the following form: 

(flu 

— = -A T Vu + V t (DVm) + Pu + /, < 7 < 1 (1) 

dP 

in which u is the scalar unknown, t denotes time, 7 is the fractional derivative 
order, V is the gradient operator, A is a coefficient vector, D is a coefficient 
matrix, P and / are scalars. Moreover, A, D, P and / are functions of 
the spatial variables. In the proposed semi-analytical FEM, the FEM is 
employed for spatial discretization which reduces the time-fractional diffusion 
equations to approximate fractional relaxation equations (also known as the 
temporal fractional ODEs). The FEM can conveniently be used to discretize 
large and geometrically complicated spatial domains. On the other hand, the 
analytical solutions for fractional relaxation equations exist and the burden 
of long-time range computation can be significantly alleviated. The present 
semi-analytical method can solve ID, 2D and 3D problems with variable 
coefficients conveniently at low implementation cost. 



2. Algorithm framework 

Time-fractional diffusion equations are often used to characterize the con- 
taminant transport processes, which exhibit typical subdiffusion features, 
such as heavy-tail decay of concentration and nonlinear time dependent mean 
square displacement (x 2 (t)) ~ t a (0 < a < 1) [4]. The time-fractional diffu- 



3 



sion problem can be expressed as: 



( (Pu 



dP 



-A T Vu + V T (DVw) + Pu + / in Q x (0, T), 

u = u on To, 

q on T N , 
h c (u - Moo) on r c , 
u\t=o = uq in Vt. 



n T (D Vw) 



(2) 



n T (DVu) 



Here, u is contaminant concentration, A is the generalized convective coef- 
ficient vector, D is the generalized diffusion coefficient matrix, Pu represents 
a reaction-absorption term, / is the source term and Q denotes the spatial 
domain of the problem. Moreover, n is the unit outward normal vector to 
the boundary, h c is the convective coefficient and r^urAfUrc = dQ which is 
denotes the entire boundary of fl The subscripts D, N and C for T designate 
the essential (Dirichlet), natural (Neumann) and convective boundary condi- 
tions, respectively, whereas r^, T N and T c are mutually exclusive. It will be 
assumed that A, D, P, f, q and independent of time. In the problem 

statement, represents the Caputo fractional derivative whose definition is 
given as below: 



<P 
dP 



9{t) 



9 r (r)di 

r(f- 7 )7 Jt^r)- 



; < 7 < 1 



(3) 



in which T is the Gamma function and 7 is the derivative order. The Caputo 



fractional derivative has the following properties |35|, |36 



P 1 

—EJXP) = XEJXP), 

f (4) 

Constant = 

dP 



in which E~ represents the Mittag-Leffler function with one parameter [20 



^) = Ef(^T)^> > *ec (5) 
The weighted residual statements for the governing equation, natural bound- 



4 



ary condition and convective boundary condition: 



ip[-j£ + A T Vu - V T (DVu) -Pu- f]dQ = 0, 



/ ^[n T (DVw) - q]dT = 0, 
Jr N 

/ ip[n T (DVu) - h c (u - u^dT = 



(6) 



can be merged to form the following weak form with the help of the divergence 
theorem: 

f ip^dtt + f ipA T (Vu)dtt + [ (V^)D(Vn)rffi - f PipudVt 
Jn dP J n Jn Jn ^ 

= / ipqdT + / iph c (u — u^dT + / ipfdVt 
Jv N Jr c Jn 

in which the trial solution u equals u and the weight function ip vanishes on 
To- In the finite element method, u and \I/ within each element fl e can be 
expressed respectively as: 



n—njy 



I'D 



u 



Ntut + R M = Nexje + N e U e 



i=l 



i=l 



n—ri£) 



(8) 



Nfifi = N e ^ e = (^ e ) T (N eV; 



where u\ and u\ are the values of u at nodes away from and on Tp, respec- 
tively; ipf are the value of %p at nodes away from Tp] Nf and Nf are the 
nodal interpolation functions for u\ and w|, respectively; n is the number of 
nodes in the element, Ud is the number of nodes on T e D = dQ e R Tp, Nfs 
and iVfs form the row interpolation matrices N e and N e , respectively; m|s 
and wfs form the vectors U e and U , respectively. By virtue of (JH]) , (JTJ) can 
be written as: 



(P_ 



or 



W([c,c]— 



u e 



u 



[K e ,K e l 



U e 

u e 



0. 



dm u 



[K,K1 



U 
U 



(9) 



(10) 
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in which 

[C e ,C e ] = / (N e f[N e ,N e ]dft, 

[K e , K e ] = / ((Nf A T [VN e , VN e ] + (VN e ) T D[VN e , VN e 



(11) 



-F(N e ) J [VN e ,VN ])(fQ- / h c (N e y [N e , N ]dT, 
F e = / q(N e ) T dT- [ h cUoo (N e ) T dT+ [ f(N e ) T dn. 

JV e N JT e c JO? 

In F e , T e N = dVt e n Y N and Y e c = dVt e n r c . Moreover, , C, C, K, K, U, 
U and F are the assembled counterparts of \I/ e , C e , C 6 , K e , K e , U e , tj e and 
F e , respectively. The arbitrary nature of \I/ leads to the following system 
equation: 



rf7 - d 1 - 

C— U + KU + C— U + KU - F = 0. 

dP dm 



If 



- d? - 

C— U + KU = F^F(t), 

dr< v 7 

the above system equation can be transformed into: 

(p ~ 

C— U + KU = 

dm 

where 

U = U + K^(KU + F - F). 
At last, the time-fractional system ([2]) leads to: 

(p ~ 

C— U + KU = 0, 

dtl 



(12) 



(13) 



(14) 



(15) 



(16) 



U 



t=o 



u . 



Exact solution of the fractional relaxation equation in (TIB"]) exists and can be 
expressed as [37|, |38[ : 



U t = EJ-MV)iJ 



(17) 



6 



where M = C X K, E 7 is the Mittag-Leffler function which has been ac 



curately evaluated by Podlubny et al |39|, |40|]. In our computations, the 
employed value of the function will be accurate up to 10~ 12 . 

It can be deduced that, if 7 = 1.0, ( ITT]) becomes the exponential solution 
of the integer order relaxation equation. Next, we decompose the Mittag- 
Leffler function in (fT7|) as: 

E 7 (-Mt 7 ) = BA t B 1 (18) 

where B is the modal matrix formed by the eigenvectors of — M. On the 
other hand, A t is a diagonal matrix whose z-th diagonal entries is E 7 (— Ajf) 
and Aj is the z-th eigenvalue of — M. Substituting (fT8"|) into ( ITTj) . we get 

U t = BA t B -1 U . (19) 

Through the above manipulations, the initial-boundary value problem in 
02]) is reduced to a initial problem through spatial finite element discretiza- 
tion. The reduced problem can be solved analytically in terms of the Mittag- 
Leffler function. The practice drastically lowers the cost associated with the 
long-time range computation of the initial-boundary value problem. 

For the essential boundary condition in f ]T3]) . the equation can be trans- 
formed as: 

C — (U - K _1 F) + K(U - K _1 F) = (20) 

which would require tj — K 1 F to satisfy a fractional relaxation equation. 
In particular, the typical case in which tj is a constant d^XJ/dt 1 = also 
satisfies ffT3l). 



3. Numerical examples 

3.1. One dimensional time-fractional diffusion equation 

The following simple time-fractional diffusion problem is considered: 

d^uixA) d 2 u(x,t) 

L = k \ x e o,L , t > 0, 

dP OX 1 (cy-i \ 

u{0,t) = u{L,t) = 0, 
u(x,0) = sin(irx/L), x G [0, L]. 
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If the diffusion coefficient k = L 2 /tt 2 , the exact solution is u exact (x,t) = 
sin(nx/L)E 1 (—t' y ) whose value at x = L/2 is shown in Figure [T] for 7 = 
0.4,0.7 and 1.0. 

To construct the spatial discretization, both linear and quadratic elements 
are attempted. For the linear element, the shape functions are: 

Ni{£) = (1 - 0/2, iv 2 (0 = (1 + 0/2, £ e [-1, 1], (22) 

and the element matrices are: 



° "6 



"21" 




1 -1 " 


1 2 




-1 1 



(23) 



in which h is the nodal spacing. For the quadratic element, the shape func- 
tions are: 



^(0 = -01 -0/2, iv 2 (o = (i-o: 

and the corresponding element matrices are: 

4 2-1 
2 16 2 
-12 4 



iv 3 (0=0i + 0A 



c- = A 

30 



3h 



7 -8 
-8 16 
1 -8 



7 



(24) 



(25) 



Table [T] lists the normalized errors at different time instants yielded by 
using 10 linear, 10 quadratic and 100 linear elements. The proposed method 
can achieve accurate results no matter linear or or quadratic elements are 
employed. As usual, the quadratic element delivers much more accurate 
results than the linear element at the same nodal spacing. Another important 
feature of this method is that the accuracy of numerical result at large time 
constants can be improved by reducing the nodal spacing h. 

To estimate the convergence ratio of the linear element and the quadratic 
element, the results in Table [2] evaluated at t — 10 but different nodal spac- 
ings are prepared and the Loo-error is: 

= max \u exact (x (i), t) -u(x(i),t)\, i = 1, 2, L/h, (26) 

It can be seen that the convergence ratio of the linear element is 0(h 2 ) and 
the quadratic element is 0(h 4 ). 

In Table [U the normalized errors increase with t. To investigate the effi- 
ciency in tackling long-time range diffusion problems, the normalized errors 
of the linear and the quadratic elements at large t are computed and listed 
in Table [3l It can be seen that the normalized errors remain fairly steady 
with respect to t. 
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Table 1: A comparison of normalized errors (Errors \(u exac t(L/2,t) — 
u(L/2,t))/u exac t(L/2,t)\) of the linear element and the quadratic element. Space 
size L = 10, time-fractional derivative order 7 = 0.8 and diffusion coefficient k = L 2 /tt 2 
in (EU). 

Time Linear element Quadratic element Linear element 



O = L/10) (> = L/10) O = L/100) 



t: 


= 





0.00000000 


0.00000000 


0.00000000 


t- 


= 


1 


1.3517e-003 


1.0525e-006 


1.3485e-005 


t: 


= 


2 


2.2855e-003 


0.4709e-006 


2.2814e-005 


t: 


= 


.3 


3.0766e-003 


1.7646e-006 


3.0726e-005 


t: 


= 


4 


3.7734e-003 


2.9057e-006 


3.7703e-005 


t: 


= 


5 


4.3983e-003 


3.9302e-006 


4.3965e-005 


t: 


= 


6 


4.9644e-003 


4.8592e-006 


4.9643e-005 


t: 


= 


7 


5.4805e-003 


5.7070e-006 


5.4825e-005 


t: 


= 


8 


5.9530e-003 


6.4838e-006 


5.9572e-005 


t- 


=0 


9 


6.3868e-003 


7.1975e-006 


6.3934e-005 



Table 2: The Loo-errors and convergence ratios of 
element (Ratio = log^oa^J L 00 ^ 2 )[log{hi/h 2 )]^ 1 
7 = 0.8, space size L = 10, diffusion coefficient k 



Nodal spacing 



the linear element and the quadratic 
41|. Time- fractional derivative order 
L 2 /ir 2 and t = 10 in fl2TJl. 



h=L/10 
h=L/20 
h=L/40 
h=L/80 
h=L/160 



Loo-error Ratio 
(Linear element) 



-error 



(Quadratic element) 



Ratio 



4.327591e-004 
1.087320e-004 
2.721688e-005 
6.806336e-006 
1.701717e-006 



7.739342e-007 

1.9928 4.909022e-008 3.9787 

1.9982 3.080541e-009 3.9942 

1.9996 1.937557e-010 3.9909 

1.9999 1.265827e-011 3.9361 



Table 3: The normalized errors (Error= \(u exac t(L/2,t) — u(L/2,t))/u exact (L/2,t)\) of 
the quadratic element. Space size L = 10, diffusion coefficient k — L 2 /ir 2 , node spacing 
h = L/100 and time-fractional derivative order 7 = 0.8 in (|21j) . 



Methods 


t=10 


t=100 


t=1000 


t=10000 


Linear element 


1.0136e-004 


8.4909e-005 


8.2652e-005 


8.2307e-005 


Quadratic element 


1.3185e-009 


1.0614e-009 


1.0242e-009 


1.0186e-009 
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3.2. One dimensional time-fractional convection- dispersion equation 

Time-fractional advection-dispersion equation (also called time-fractional 
Fokker-Planck equation), which exhibits heavy-tail concentration decay fea- 
ture, is usually used to characterize contaminant transport in natural porous 
media. A simple example is: 



(Pu(x,t) du(x,t) , d 2 u(x,t) 

-^r 1 = - a ^T + k —5x^' ^(o,l), t>o, 

u(0, t) = E 7 (-(a - k)t r ), u{L, t) = e L E 1 {-{a - k)t~<), 
u(x,0) = e x , x G [0,L]. 



(27) 



Assuming a and k are constants, a > k, the exact solution of above equation 
can be written as: 



u e xact(x,t) = e x E 1 (-(a - k)^) 



(28) 



which is portrayed in Figure [2] for t < 1. For the linear element, the corre- 
sponding element matrices are: 



h 
6 



2 1 
1 2 



1 1 
•1 1 



k 

+ h 



1 -1 
-1 1 



(29) 



The normalized errors obtained by using 10, 20 and 40 elements at x — 
L/2 and different t are shown in Table HI With only 10 elements, the errors 
have been less than 0.1%. 



3.3. Two dimensional time-fractional diffusion equation 

In this example, the following two-dimensional problem is considered: 

(30) 



d?u{x, y, t) 
dB 



k(- 



dx 2 dy 2 
u(x, y, t) = 0, (x, y) G dil, t > 0, 
u(x, y, 0) = sin(xTr / L)sin(yn / L) , (x, y) G Q U dQ, 



in which k is the diffusion coefficient, Q — [0, L] x [0, L\. If k = 1/ir 2 and L = 
1.0, the exact solution of the problem is u exact (x, y, t) = sin(x7r)sin(yir)E 7 (— 2t 7 ). 

The square problem domain is modeled by 4 x 4, 8 x 8 and 16 x 16 
four-node square elements. The element interpolation functions are: 



N t = (1 - 0(1 " W4, N 2 = (1 + 0(1 - r/)/4, 



l + 0(l + 77)/4, Ar 4 = (l-0(1 + W4, S,V£[-1, 1] 



(31) 
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Table 4: The normalized errors (Error= \(u exact (L/2, t) ~ u(L/2, t))/u exac t(L/2, t)\) of the 
linear element. Space size L = 1.0, diffusion coefficient k = 1.0, convective coefficient 
a = 2.0 and time-fractional derivative order 7 = 0.8 in (|27|) . 



Nodal spacing 


t=2.0 


t=4.0 


t=6.0 


t=8.0 


h=L/10 


0.9860e-004 


0.9790e-004 


0.9724e-004 


0.9677e-004 


h=L/20 


0.2459e-004 


0.2441e-004 


0.2425e-004 


0.2413e-004 


h=L/40 


0.6143e-005 


0.6099e-005 


0.6058e-005 


0.6029e-005 



Table 5: The normalized errors (Error= \(u exac t(L/2, L/2, t) — 

u(L/2, L/2, t))/u exact (L/2, L/2, t)\) of the four-node square element. Space size 
L = 1.0, nodal spacings h = h x = h y , diffusion coefficient k = 1/tt 2 and time- fractional 
derivative order 7 = 0.8 in (|5D)) . 



Space step 


t=2.0 


t=4.0 


t=6.0 


t=8.0 


h=L/4 


6.5673e-002 


6.1143e-002 


5.7924e-002 


5.6111e-002 


h=L/8 


1.6937e-002 


1.5786e-002 


1.4929e-002 


1.4444e-002 


h=L/16 


4.2664e-003 


3.9778e-003 


3.7602e-003 


3.6368e-003 



and 



36 



4 2 12 
2 4 2 1 
12 4 2 
2 12 4 



k 
6 



4 


-1 


-2 


-1 


-1 


4 


-1 


-2 


-2 


-1 


4 


-1 


-1 


-2 


-1 


4 



(32) 



Following the calculation steps ( jl6]) - (jl~9l) . Figure [3] plots the numerical 
solution for 7 = 0.8att = 2, the normalized errors for 7 = 0.8 at x = 
L/2, y = L/2 and various values of t are given in Table The errors 
drop with the nodal spacings. Indeed, the finite element method can readily 
take coordinate-dependent and direction-dependent diffusion coefficients into 
account. 



3.4- Time-fractional diffusion equation in a quarter of circular domain 

An important advantage of the finite element method over the finite dif- 
ference method is that the former can readily consider complex spatial do- 
mains. In this example, the following problem defined over a circular domain 
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is considered: 



<Fu{x,y,t) _ u( d 2 u{x,y,t) d 2 u(x,y,t) ^ 

n T (DV«) = 0, x = 0,y^l, and y = 0, < x < 1, t > 0, (33) 
u(x,3/,t) = J (l)£ 7 (-t 7 ), G {(x,y)|x 2 + y 2 = 1}, 



0) = J (^ 2 + 2/ 2 ), 0, y)enudtt 



in which fi = {(x,y)\x > 0, y > 0, x 2 + y 2 < 1}. If k = 1, the exact solution 
of ( )33l) is u exact (x,y,t) = Jo{\/ x 2 + y 2 )E^(— t 7 ), in which J represents the 
zeroth order Bessel function of the first kind. For symmetry, we only need 
to model a quarter of the problem domain and a typical mesh is depicted in 
Figure HI 

For the four-node element, the (x,y) coordinates are also interpolated 
with the functions given in (l3Tj) . i.e. 



i=i i=i 
in which yf) are the coordinates of the i-th element nodes. 




i ,i 

/ (N e ) T N e (iet(J)ci^, 
-l J-i 

(VN e ) T J- 1 DJ- T (VN e )rfet(J)d^- 

u-i 



(35) 



d/di 



and J 



dx e /d^ dy e /d£ t 
dx e /dr) dy e /dr] 



(36) 



The matrices C e and K e are computed by the second order Gauss- Legendre 
rule. In our computations, a quarter circle is partitioned into 3, 48 and 217 
elements, the corresponding numerical results at some selected spatial loca- 
tions are listed in Table [6] for 7 = 0.8. Table [6] indicates that the proposed 
method is capable of delivering accurate solution to anomalous diffusion prob- 
lem ( ]33|) and the accuracy can be improved by employing more elements in 
modeling the computational domain. 
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Table 6: The numerical results by different numbers of elements, In the computation, 
k(x) = k(y) = 1, the exact solution of (130) is u exact (x,y,t) = J (\/x 2 + y 2 )E-y(~ i 7 ), 7 = 
0.8. 



Coordinates 


3 Elements 


48 Elements 


217 Elements 


Exact solution 


(0,0) 


0.37770 


0.38638 


0.38681 


0.38695 


(0.35355, 0.35355) 


0.34055 


0.36233 


0.36292 


0.36314 


(0.21339, 0.21339) 




0.37760 


0.37804 


0.37819 


(0.42678, 0.17678) 




0.36603 


0.36644 


0.36658 


(0.67533, 0.27973) 




0.33659 


0.33687 


0.33696 


(0.53033, 0.53033) 




0.33404 


0.33432 


0.33442 


(0.27973, 0.67533) 




0.33659 


0.33687 


0.33696 


(0.17678, 0.42678) 




0.36603 


0.36644 


0.36658 



4. Application 

To investigate the efficiency and applicability of the proposed semi-analytical 
method in solving real-world problems, it is employed to solve the problem of 
tracer solute transport in an aquifer. The experiment was conducted using a 



test aquifer in Nevada and the schematic diagram is shown in Figure [5] 13 
Bromide of quantity M = 20.81 kg was used as a nonsorbing tracer solute 
and introduced to the injection well for a period of T = 3.54 days at an 
average concentration of 3.77 kg/m 3 . A reference point, the injection well 
and extraction well are located at r = m , 30 m (Ri) and 60 m (R e ) along 
the downstream direction of the underground water flow. The radius of the 
extraction well is 0.127 m, the center of extraction well is at r c = 60.127 m, 
the pumping rate is Q = 12 A m 3 /d. The solute concentration in the extrac- 
tion well had been monitored for about 321 days and the screened interval 
was b = 35 m . More detailed description of the experiment can be found in 
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Since To is short compared with the total time range (~ 321 days) of 
measurement, the following radial initial-boundary value problem for the 
solute transport in the fractured granite aquifer is established in terms of the 
solute concentration u as: 

( (Pu{r,t) v du(r,t) 1 d du(r,t) 

1 Tr{ a o — 7; J; r t \y,rL e ), 



dP Tc. — t dr r c — r dr dr 

du(Re,t) = f 

or 

u(r,0)=/(r), re [0, BJ 



= 0, *gJ)=o, t >o, < 37) 
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where v is the convective coefficient, d is the dispersion coefficient, v = ad 
and a is the dispersivity. Moreover vq/ (r c —r) and do/(r c —r) have the units of 
[L/T 7 ] and [L 2 /T 7 ] which represent the nonlocal aquifer properties jjjj. The 
initial value is normalized as f(r) = M5(r — (r c — Ri))/(2ir(r c — Ri)b9T ), 6 
is the hydraulic parameter. The boundary conditions imply that the solute 
cannot reach r = by upstream dispersion and the solute moves by advection 
at r = R e which gives the wall of the extraction well (l3| . 

In order to obtain a high accurate numerical approximation, the quadratic 
element is adopted, the element matrices C e and K e are computed by the 
second order Gauss-Legendre rule. 

A comparison of the numerical predictions and the experimental data 
is shown in Figure and the heavy-tail feature characterized by the time- 
fractional model ( 157)) with different derivative orders is shown in Figure [7J 
It can be observed from Figure |6] that the numerical result offers a good fit 
to most of the experimental data. Due to the subdiffusion behavior in the 
aquifer matrix and immobile water, in the experimental result, the concen- 
tration of bromide exhibits a rather slow decay in the late time. Figure [7J 
confirms that the time-fractional radial flow model ( 157)) captures the long- 
time behavior with heavy-tail. Figure UJ also illustrates that the heavy-tail 
feature becomes more remarkable with the decreasing of the time-fractional 
derivative order 7. Hence, in this model (157)) . the time- fractional derivative 
order 7 is a indicator of the non-Fickian transport caused by the complex 
structure of the fractured aquifer. 

5. Discussions 

By using the finite element to discretize the spatial domain, the fractional 
diffusion equations can be reduced to approximate fractional relaxation equa- 
tions which possess analytical solutions. The semi-analytical method can not 
only compute time-fractional diffusion equations in long-time range at low 
computational cost but also deliver accurate numerical predictions for com- 
plex and large spatial problem domains. The accuracy in spatial domain can 
be improved by using more elements, high-order elements or elements based 
on advanced finite element formulations. Since the exact solution is used 
in time domain, the stability and convergence conditions of the proposed 
method can be easily satisfied. It can be said that the proposed method is 
more robust than previous ones. 
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The main restriction for the proposed method is that the weak forms of 
time diffusion equations can be transformed into the following form: 

C^u + Ku = 0. (38) 

In cases that the source term, physical parameters and/or boundary condi- 
tions are only weak function(s) of time, a multiple time step method can be 
used. 



6. Concluding remarks 

From formulations and examples presented, it is clear that a class of 
time-fractional diffusion equations can be easily computed and the heavy- 
tail feature can be accurately characterized by the new method. Our future 
research work will focus on advanced finite element formulations, such as hy- 



brid element [461 ]. to compute temporal-spatial fractional diffusion equations 



which characterize more complex contaminant transport problems. 
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Figure 1: The diffusion curves of time-fractionai diffusion modei (f2Tj) at x = L/2, obtained 
by the quadratic eiement. Nodai spacing h = L/100 and space size L = 10. 
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a=2.0, k=1 .0, exact solution 
+ a=2.0, k=1 .0, numerical solution 

— a=4.0, k=1 .0, exact solution 

o a=4.0, k=1 .0, numerical solution 

— a=2.0, k=1 .5, exact solution 

❖ a=2.0, k=1 .5, numerical solution 



Figure 2: A comparison of exact and numerical solutions of time- fractional advection- 
dispersion model (|27|) at x = L/2 with 7 = 0.8. Nodal spacing h — L/20 and space size 
L = 1.0 in the linear element. 
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Figure 3: The numerical result of two dimensional time-fractional diffusion model pop 
with 7 = 0.8 at t = 2.0. Diffusion coefficient k = 1/ir 2 , space sizes L = 1.0 and nodal 
spacing h = L/16. 
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Figure 4: The quadrant of a circle domain with unit radius meshed into 48 elements. 
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Figure 5: A schematic diagram of experiment. 
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Figure 6: The numerical approximation of the time- fractional radial flow advection- 
dispersion equation model (|33[) with 7 = 0.92 at r = R e . In this numerical simulation, con- 
vective coefficient vq = 0.0564/0, 9 — 0.023, the dispersion coefficient do = a^c a = 6.8, 
nodal spacing h = 3 and time step is At = 10. 
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Figure 7: The numerical approximation of the time-fractional radial flow advection- 
dispersion equation model Q33[) with different time- fractional derivative values 7 at r = R e . 
In this numerical simulation, convective coefficient vq = 0.0564/6*, 9 = 0.023, dispersion 
coefficient do = avo, a = 6.8, nodal spacing h = 3 and time step is At = 10. The curve 
with 7 = 1.0 corresponding to Fickian dispersion. 
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